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1. Introduction 

a number of models have been published that relate different physiological processes in¬ 
volving glial cells to calcium dynamics. De Pitta et al. [22] give an overview of current prob¬ 
lems in the modelling of astrocytes. One area of continuing interest is the propagation of 
signals between astrocytes via intercellular calcium waves. Hofer et al. [42] investigated the 
spreading of signals between astrocytes via calcium waves based on a model by Sneyd et al. 
[76]. Bennett et al. [8, 7] developed a more detailed model of calcium waves that combines un¬ 
derlying calcium dynamics with ATP release by purrnergic receptors in order to demonstrate 
that calcium waves depend on ATP release rather than on IP3 diffusion through gap junctions 
as in the model by Hofer et al. [42]. Edwards and Gibson [27] later published a model that in¬ 
cluded both modes of signal propagation and concluded that both were necessary to account 
for data collected from the retina. Recently, the study of calcium waves has been extended from 
one- or two-dimensional to three-dimensional spatial domains [45]. Macdonald and Silva [50] 
model wave propagation on an astrocyte network derived from experimental data. The Ben¬ 
nett et al. model was used for inveshgahng spreading depression, a wave of electrical silence 
that propagates through the cortex and depolarises neurons and glial cells [9]. 

A fundamental problem in calcium dynamics in general is the question how multiple signals 
can be encoded by the dynamics of a single quantity, the concentration of calcium. De Pitta 
et al. [24, 20, 23] investigated how a stimulus could be encoded via the frequency or the am¬ 
plitude or both frequency and amplitude which demonstrates that two different signals can 
be represented independently in an individual calcium signal. Dupont et al. [26] showed in 
a detailed model how the signal received by a particular glutamate receptor is encoded via 
calcium oscihahons. 

Lavrentovich and Hemkin [46], Zeng et al. [87], Riera et al. [62,63] inveshgated spontaneous 
calcium oscillations in astrocytes and Li et al. [47] explored their role in spreading depression. 

Also the coupling of astrocyte network with the neural network has been investigated. At 
the single-cell level, De Pitta et al. [21] modelled the interaction of an astrocyte with a synapse. 
Allegrini et al. [1], Postnov et al. [59] study the influence of a network of astrocytes on a neural 
network. 

Most recently. Barrack et al. [5, 6] explored the role of calcium signalling in neural develop¬ 
ment. By coupling calcium dynamics with a model of the cell cycle they examine how glial 
progenitors differentiate to neurons triggered by a calcium signal. 

This review of the modelling literature on glial cells clearly demonstrates that the impor¬ 
tance of calcium d}mamics is well recognised—the majority of studies in the literature accounts 
for calcium signalling and often models are used to find a link of physiological processes with 
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calcium signalling. In many cell types including glial cells the inositol trisphosphate receptor 
(IP 3 R) plays a crucial role in inducing oscillatory Ca^^ signals. In the presence of IP 3 , opening 
of IP 3 R channels leads to Ca^+ release from the endoplasmic reticulum (ER), an intracellu¬ 
lar compartment with a very high Ca^+ concentration a few orders of magnitude higher than 
that of the cytoplasm. The IP 3 R is activated by Ca^+ so that such a release event dramatically 
increases the open probability of the IP 3 R which induces further release of Ca^+(henceforth 
called calcium-induced-calcium-release, or CICR) until a high Ca^^ concentration in the chan¬ 
nel environment eventually inhibits the IP 3 R. 

The Li-Rinzel model [48], an approximation of the classical De Young-Keizer model [25], is 
by far the most commonly used representation of the IP 3 R in models of glial cells. Only Alle- 
grrni et al. [1] and Lavrentovich and Hemkrn [46] chose different models based on Atri et al. 
[3] or Tu et al. [82], respectively. Dupont et al. [26] use the model by Swillens et al. [79] that 
explicitly accounts for the effect of interactions in a cluster of IP 3 R channels. Early models of 
the IP 3 R were designed to account for the bell-shaped Ca^+ dependency of the open probabil¬ 
ity po of the channel described by Bezprozvanny et al. [10]. Since then the dynamics of IP 3 R 
in response to varying concentrations of IP 3 , Ca^^ and ATP has been characterised much more 
comprehensively as well as the differences between the different isoforms of the IP 3 R (among 
the models mentioned above, in fact, only Tu et al. [82] accounts for the fact that astrocytes 
predominantly express type IIIP 3 R). 

The scope of current data-driven models of ion channels has advanced beyond representing 
the average open probability po- Recent models capture the stochastic opening or closing of 
single IP 3 RS in aggregated Markov models i.e. instead of only modelling the stationary be¬ 
haviour of the channel they represent the dynamics of the IP 3 R (Section 3.4). Accurate repre¬ 
sentation of IP 3 R dynamics depends on various sources of experimental data (Sections 3. 1-3.2) 
as well as appropriate statistical methods for fitting Markov models to these data (Section 3.5). 
Statistical methods automate the process of estimating parameters for a given Markov model. 
Thus, the main challenge of data-driven ion channel modelling is to define the structure of 
a Markov model which allows the integration of various sources of experimental data. We 
illustrate this process with two recent examples of models for the IP 3 R (Sections 3.6 and 3.7). 

Once a model for a single channel has been developed, data from small clusters of channels 
can be used to determine how well the behaviour of a cluster is represented by an ensemble 
of single-channel models (Section 4.1). Studying the influence of an IP 3 R model on calcium 
dynamics allows us to evaluate the relative importance of different aspects of single-channel 
dynamics. Cao et al. [14] showed that the essential features of calcium dynamics in airway 
smooth muscle could be preserved after iteratively simplifying the IP 3 R model by Siekmann 
et al. [73] to a two-state model that only accounted for the switching between the inactive 
"park" and the active "drive" mode. In Section 4.2 it is shown that this also applies to the puff 
distribution. This demonstrates that modal gating is the most important regulatory mecha¬ 
nism of the IP 3 R. It also emphasises that data-driven modelling of ion channels does not nec¬ 
essarily have to lead to detailed models based on complicated model structures but rather can 
be used so that relevant data is selected to represent ion channels at the appropriate level of 
complexity for a given application. 


2. Mathematical models of calcium dynamics/ CICR 

The purpose of a mathematical model of CICR is to explain the emergence of complex intra¬ 
cellular calcium dynamics such as oscillations as the result of interdependent calcium fluxes. 
This comprises both fluxes into and out of the cell as well as the exchange between the cy¬ 
tosol and intracellular stores (Pigure 1). The dynamics of cytosolic (c) and stored calcium (cer) 
resulting from these fluxes can be represented by a system of differential equations: 
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Figure 1. General structure of calcium fluxes in glial (and other non-excitable) 
cells. The central component is the flux /ipr through the inositol trisphosphate 
receptor (IP 3 R). The IP 3 R is activated by binding of IP 3 which is generated upon 
stimulation of the cell by an agonist. This causes the release of Ca^+ from the 
endoplasmic reticulum (PR) to the cytoplasm. The resulting elevated Ca^^ con¬ 
centration increases the open probability of the IP 3 R and the ryanodine recep¬ 
tor (RyR) which stimulates further Ca^^ release. This mechanism is known as 
calcium induced calcium release (CICR). At high concentrations, Ca^+ inhibits 
the IP 3 R, i.e. the open probability of the IP 3 R decreases. In consequence, /serca 
influx into the PR through the SERCA pump dominates the efflux through IP 3 R 
and RyR so that Ca^^ is reabsorbed by the PR. Ca^+ exchange with the extra¬ 
cellular space is controlled by uptake through various channels (/in) and by 
extrusion via pumps (/pm)- 


( 1 ) 

( 2 ) 


dc 

It 

rfCER 

dt 


/1P3R + /RyR + Jin — /pm — /sERCA 

7 (/serca — /1P3R — /ryr) 


Here, /in is calcium influx from the extracellular space via calcium channels located in the cell 
membrane, /pm accounts for calcium removed from the cell by the plasma membrane pump. 
/ 1 P 3 R and /RyR represent calcium release from the endoplasmic reticulum (PR) through the IP 3 R 
and the RyR, respectively, and /serca stands for reuptake of calcium into the PR by the SERCA 
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pump. The conversion factor 7, the ratio of the cytoplasmic volume to the ER volume, is nec¬ 
essary because calcium concentrations are calculated with respect to the different volumes of 
these two compartments. The model (1), (2) provides a description of Ca^+ concentrations 
across the whole cell. This means that we cannot account for spatial effects due to hetero¬ 
geneities of the spatial distribution of IP 3 R, SERCA and other relevant components of the sys¬ 
tem. By using a deterministic model we further assume that the various Ca^+ fluxes can be 
described as deterministic after averaging over a large number of channels and transporters. 
In Section 4 we will consider a stochastic model over a small spatial domain for a cluster of 
interacting IP 3 RS. 

In a whole-cell model of calcium dynamics such as (1), (2), a representation of the IP 3 R must, 
in principle, just provide a functional expression for 


(3) /iP3R([IP3],[Ca2+],[ATP]), 

the ligand-dependent flux through IP 3 R channels present in a cell. Because the calcium con¬ 
centration [Ca^+] is time-dependent, /ipjR varies over time. In the early days of modelling of 
the IP 3 R, phenomenological models were used for representing the IP 3 R flux. A good example 
is the model by Atri et al. [3]: 


(4) 7,p,r(p,c) = N„pe„t (w + 

where p = [IP 3 ], c = [Ca^+] and Nopen is the number of open channels. The model by 
De Young and Keizer [25] is derived from more detailed assumptions on chemical interactions 
of the channel with its ligands. In Section 3.6 we present a more recent model [84] that is 
representative for this approach. The Hill function-type terms in (4) enabled Atri et al. to 
interpret their model in terms of a physical process but the main motivation of the model was 
to obtain a fit of the calcium-dependent whole-cell flux /iPgR to data collected by Parys et al. 
[58]. Erom a purely mathematical point of view, phenomenological models seem to be the 
ideal approach for investigating the role of IP 3 R in calcium dynamics—restriction to minimal 
models that generate the desired behaviour ensures that model behaviour can be analysed to 
a great extent. This allows us to test hypotheses on IP 3 R regulation in an elegant way. 

But the capability of simple mathematical expressions for the macroscopic flux /iPgR to per¬ 
form the appropriate functional role in calcium dynamics is only a relatively indirect test for 
IP 3 R models. By following a phenomenological approach we mostly ignore data that gives 
more direct information on the IP 3 R, such as, for example, the molecular structure of the chan¬ 
nel protein which can be obtained from crystallography and time series of opening and closing 
of a single channel from patch-clamp recordings. Taking into account these data may allow us 
to restrict the set of theoretically possible mathematical expressions and, in this way, also the 
set of possible mechanism. 

3. Data-driven modelling of single IP 3 RS 

Because most biophysical data relate to single channels, data-driven modelling involves an 
important conceptual step—instead of directly specifying the whole-cell flux / 1 P 3 R, we first 
construct a model for the flux through a single channel. Whereas for the macroscopic flux /iPgR 
which is averaged spatially over many channels distributed across the whole cell the deter¬ 
ministic model (3) is appropriate, representing the flux through a single channel requires a 
stochastic model. In a second step, /iPgR is then derived by appropriately averaging over the 
stochastic fluxes through individual channels. 
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In Sections 3.1 and 3.2 we describe two sources of data that are commonly used for the con¬ 
struction of ion channel models. Ca^^ release data from small clusters of IP 3 R, so-called cal¬ 
cium puffs (Section 3.3), can be used for validating models of single channels. In Section 3.4 ag¬ 
gregated continuous-time Markov models, the mathematical framework common to all mod¬ 
els based on single-channel data, is introduced. A short review of statistical approaches for 
fitting Markov models to single-channel data is given in Section 3.5. In Sections 3.6 and 3.7 
examples of two recent models of the IP 3 R are given in order to illustrate different modelling 
approaches. Earlier models have been reviewed by [38] and Sneyd and Falcke [77]. Model 
comparisons [78, 40] generally show that models not parameterised by fitting to experimen¬ 
tal data may not do a very good job at reproducing the statistical properties of ion channel 
kinetics. 


3.1. Molecular structure. The mathematical structure of many ion channel models is designed 
to mimic the chemical structure of the channel protein. The motivation for this approach is to 
link molecular structure of the ion channel to its function. 

In vertebrates there exist three different genes encoding three different types of the IP 3 R. In 
mammals, type IIP 3 R is ubiquitously expressed but most cells express more than one isoform. 
The predominant isoform in astrocytes is type II IP 3 R [69, 43]. For each isoform there are 
several splice variants. 

Imaging the three-dimensional structure of the complete IP 3 R and RyR channel proteins is 
challenging and only recently have accurate 3D visualisations of complete IP 3 RS using electron 
cryomicroscopy (cryo-EM) become available [49]. Parts of the channel can be imaged at higher 
resolution by crystallography and be superimposed on cryo-EM images [29]. These studies 
have revealed that IP 3 R channels are tetramers i.e. formed by binding of four IP 3 R proteins. 
These tetramers may consist of different IP 3 R subtypes but experimental studies have so-far 
concentrated on investigating homotetramers formed by four copies of the same subtype (but 
see Alzayady et al. [2]). The classical model by De Young and Keizer [25] took into account 
this information by building a model from identical subunits that all had to be in an open state 
for the channel to open (although the de Young-Keizer model assumed three instead of four 
subunits). 

Analysis of the amino acid sequence by mutation experiments have assigned functional 
roles to various segments, for example, the IP 3 binding core (IBC) which contains an IP 3 bind¬ 
ing site has been identified. There is less information on the number and localisation of Ca^^ 
binding sites. Because localisation of Ca^+ binding sites by mutation studies has been dif¬ 
ficult, Foskett et al. [31] infer various Ca^+ binding sensors from the observed co-regulation 
by IP 3 and Ca^+, see Foskett and Mak [30] for a summary. Often models assume a certain 
number of IP 3 and Ca^+ binding sites and represent binding and unbinding of these ligands 
as transitions between states regulated by mass action kinetics. This modelling approach will 
be described in more detail in Section 3.6. 


3.2. Patch-clamp recordings. Detailed studies of individual ion channels became possible due 
to the development of the patch-clamp technique. Neher and Sakmann [56] were the first to 
detect the flow of ions through a single ion channel by measuring the resulting current at 
constant voltage. The time-course of opening and closing can be inferred from the detected 
current which stochastically jumps between zero (closed) and one or more small non-zero 
current levels in the pA range (open) whose sign depends on the valence of the ion and the 
direction of the current. 

Mak and Foskett [53] recently reviewed the single-channel literature of IP 3 R channels. An 
important experimental development that they highlight relates to the difficulty that IP 3 RS are 
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naturally localised within cells rather than in the cell membrane. Whereas in earlier patch- 
clamp experiments, IP 3 R channels were studied in artificial lipid bilayers, more recently in¬ 
vestigating IP 3 R in isolated nuclei is favoured because it is assumed that nuclei provide an 
environment similar to the endoplasmic reticulum (ER), the native domain of the IP 3 R. 

3.2.1. Stationary data. If ligand concentrations (such as IP 3 , Ca^+ and ATP) are kept constant 
for the whole duration of the experiment we obtain stationary data. These data allow us to 
observe the "typical" channel dynamics for a given combination of ligands. The reason that 
we refer to these data as "stationary" is that we assume that the channel has fully adjusted to 
the concentration of ligands—the term stationary suggests that the channel has reached its sta¬ 
tionary probability distribution, see Section 3.4. Because the stationary solution is only reached 
asymptotically we can, in theory, never be sure that our ion channel has actually reached equi¬ 
librium. Instead we can check if a data set is not stationary by using indicators such as the open 
probability. If the open probability averaged over a sufficient number of data points sponta¬ 
neously changes (which indicates the switching of the channel to a different activity level) the 
channel may exhibit modal gating. 

3.2.2. Modal gating. Spontaneous switching between different levels of channel activity at con¬ 
stant ligand concentrations has been observed for a long time. The earliest example is perhaps 
from a classical study of the large-conductance potassium channel (BK) [52, 51]. In IP 3 R chan¬ 
nels modal gating was discovered only relatively recently [44]. The authors found three differ¬ 
ent modes characterised by high (H), intermediate (I) and low (L) levels of open probabilities, 
Pq, Pq and Pq. They also realised the importance of modal gating for IP 3 R regulation: they 
observed that the same three modes seemed to exist for different combinations of ligand con¬ 
centrations. Because the IP 3 R mostly seemed to adjust the time spent in each of the three 
modes they proposed that modal gating is the major mechanism of ligand regulation in IP 3 R 
channels. 

One reason that the significance of modal gating has not been appreciated until recently 
is due to the fact that switching between different modes cannot always be recognised eas¬ 
ily without statistical analysis. Recently, Siekmann et al. [72] developed a statistical method 
which for a given set of single-channel data detects switching between an arbitrary number 
of modes M' characterised by their respective open probabilities Pq'. A software implemen¬ 
tation which is publicly available under https: //github. com/merlinthemagician/icmcstat. 
git was applied to a large data set from Wagner and Yule [ 86 ]. Siekmann et al. [72] found that 
the same two modes, an inactive "park" {p^^^ ~ 0 %) and an active "drive" mode (Pq~ 
70%), were found across all combinations of ligands. There may be various reasons why two 
modes were observed rather than the three modes found in the earlier study [44], see Siek¬ 
mann et al. [72] for more details. But more importantly, a detailed study of a bacterial potas¬ 
sium channel (KscA) [17, 16, 15] strongly suggests that the stochastic dynamics characteristic 
for each mode may be closely related to distinct three-dimensional configurations (conforma¬ 
tions) of the channel. Thus, whereas it is often difficult to relate individual open or closed 
states in ion channel models to distinct conformations of the channel protein, the set of model 
states that represents a particular mode may, in fact, have a biophysical counterpart [72]. In 
order to confirm this hypothesis, more studies of modal gating for a variety of channels are 
needed. 

Independent from its biophysical significance, appropriately accounting for modal gating is 
crucial from a modelling point of view. As we will see in Section 3.4, the phenomenon of modal 
gating demonstrates that a Markov process must be observed for a sufficiently long time in 
order to infer the correct stationary distribution, otherwise we observe a "quasi-steady state". 
For example, a channel whose kinetics is restricted to an active and an inactive mode can 
produce intermediate activity only by switching between both modes. Thus, a model that is 
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not capable of switching between different levels of activity is misleading because it produces 
a constant open probability instead of alternating between highly different open probabilities. 
In their recent review Mak and Foskett [53] explicitly recognise the importance of modal gating 
which so far has only been represented in the most recent models [84, 73]. 

3.2.3. Response to rapid changes of ligand concentrations. Modal gating is an aspect of stationary 
data collected at constant concentrations of ligands. In contrast, Mak et al. [54] designed an 
experiment where IP 3 and/or Ca^+ concentrations in the channel environment were rapidly 
altered in order to simulate an instantaneous change of ligand concentrations. Switching from 
ligand concentration where the IP 3 R is inactive to conditions where the channel is maximally 
activated (and vice versa), enabled Mak et al. [54] to investigate the question how fast the IP 3 R 
responds to varying ligand concentrations. To illustrate the experiment let us consider the 
change from inhibitory to activating conditions. At an inhibitory condition, the open prob¬ 
ability of the channel is very close to zero (po ~ 0) so that initially the IP 3 R is most likely 
closed. When changing from an inhibitory to an activating condition the channel will activate 
but it needs a certain time to respond to the change. In order to measure this latency, Mak 
et al. [54] recorded the time the channel took from when they altered the ligand concentration 
until the first opening. For the opposite change from activating to inhibitory conditions they 
analogously detected the time the channel needed to switch from a high to a low level of activ¬ 
ity. This experiment was repeated multiple times for switching between the same conditions 
which enabled the authors to investigate the latency statistics. It was not only discovered that 
for some conditions the latencies were surprisingly long but interestingly, they also found that 
for some conditions the latency distributions were multi-modal which shows that multiple 
timescales may be observed for the same latency. 

Due to the substantial effort required to perform these experiments (which have to be re¬ 
peated multiple times for each condition where each repeat only gives a single data point 
rather than a time course) it is unsurprising that these data are very rare. In fact, to date, Mak 
et al. [54] is the only data set of this kind currently available. Mak and Foskett [53] explain that 
their data suggests that there may be long refractory periods between release events from the 
same IP 3 R which makes these results particularly relevant for the modelling of Ca^^ puffs. 

3.3. Calcium puffs. So far we have only considered data recorded from single IP 3 RS. In or¬ 
der to understand how the macroscopic flux / 1 P 3 R arises from the release of many individual 
channels we have to consider the hierarchical nature of Ca^+ release. As reviewed by Parker 
et al. [57], Falcke [28], Thurley et al. [80] stochastic opening of a single IP 3 R channel leads to 
a localised Ca^^ release event (a Ca^^ blip). Such a release further sensitises neighbouring 
IP 3 R to induce more Ca^+ release through a few tightly clustered IP 3 R by CICR (a Ca^+ puff). 
Sufficiently many puffs could eventually trigger a global elevation of [Ca^+]i that is able to 
propagate through the entire cell (a Ca^+ wave) [55]. Thus, Ca^+ puffs play a crucial role that: 
not only are puffs essential for the formation of functional global Ca^+ signals [ 12 ] but they 
also reflect the quantal Ca^+ releases by stochastic openings of IP 3 R in vivo [75]. 

Experimentally, Ca^+ release at a specific spatial position can be initiated by triggering re¬ 
lease of caged IP 3 using a laser. A relative measure for the local Ca^+ concentration is obtained 
by detecting fluorescent dye bound to Ca^^ using a light microscope. For a given point within 
the cell the resulting time series is characterised by a sequence of stochastic spikes that are 
highly variable as far as the spike amplitude, the frequency and the time interval between sub¬ 
sequent spikes, the inter-puff interval, is concerned. From a modelling point of view, these 
data can be used to test wether the single-channel behaviour represented in a model is able 
to account for the release from a cluster of interacting IP 3 RS. As explained in Section 4.1, Cao 
et al. [13] found that the original model by Siekmann et al. [73] was incapable of generating the 
correct stochastic puff distribution as long as the adaptation to different ligand concentrations 
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was assumed to occur instantaneously. After augmenting the model so that it accounted for 
the latency data by Mak et al. [54] presented in the previous section the puff statistics could be 
reproduced accurately 

The only other model that accounts for latency data is the model by Ullah et al. [84]. Because 
the models by Siekmann et al. [73], Cao et al. [13] and by Ullah et al. [84] are the only models 
that account for all aspects of single-channel data assumed to be necessary for an understand¬ 
ing of the IP 3 R we focus on these two models and the alternative modelling approaches that 
they represent in Sections 3.6 and 3.7. 

3.4. Aggregated continuous-time Markov models. The most natural model for the stochas¬ 
tic process of opening and closing of a single ion channel is the aggregated continuous-time 
Markov model. A good introduction to the theory reviewed here is the classical paper by 
Colquhoun and Hawkes [18] which also gives some simple but illustrative examples. 

An aggregated continuous-time Markov model is a graph on a set of tiq closed and no open 
states S = {Cl,... ,Cnc,0„^+i,... ,Onc+no} (Figure 2). 




(a) Ullah et al. [84] (b) Siekmann et al. [73] 

Figure 2. Examples for recent Markov models of the IP 3 R. 

Between adjacent states S; and Sj the transition rate (from S; to Sj) is given by qij > 0 so 
that the whole model is represented by a matrix with constant coefficients, the infinitesimal 
generator Q = [qij). The time-dependent probability distribution p(t) over the state set S is 
the solution of the differential equation 

(5) = p{i)Q, P(0) = po- 

The stochastic interpretation of (5) is as follows: for a given point in time, one particular 
state Si of the model is "active". But how long it will take until the current state S, is vacated 
and which state Sj will be active after a time t cannot be answered with certainty (determinis¬ 
tically) due to the stochastic transitions between states. 

For the model defined by (5) the Markov property holds both for the stochastic sequence of 
active states as well as for the time that it takes until the active state is left. 
























(1) Which state Sj will be the next active state only depends on the currently active state Si, 
not on previously active states. 

(2) The time ts, h takes until the model exits from the state Si, also called the sojourn time 
in Si, does not depend on the time already spent in S,. 

The second point implies that sojourn times fs, must be exponentially-distributed because 
the exponential distribution is the only continuous probability distribution with this property 
This explains why multiple open and closed states may be needed for accurately representing 
the opening and closing of ion channels. 

In order to ensure that p{t) is a stochastic vector i.e. Yllh Vi’ for f > 0/ the matrix Q 

must be conservative, i.e. for the diagonal elements qu we have 

(6) qa = - X] ‘ly' i,; = 1,..., ng. 

Provided that (6) holds, the solution 


(7) p{t) = poexp{Qt), 

is a stochastic vector for all f > 0 if and only if the initial distribution po is a stochastic 
vector. From (7) the time-dependent open probability po(0 of the charmel can be calculated 
by summing over the individual probabilities of all open states. 

For large times t the solution p(f) approaches a stochastic vector n which is known as the 
stationary distribution. This means that provided we wait sufficiently long, the expected fre¬ 
quency of observing a state Si approaches a probability tt;. Because p{t) is the solution of a 
differential equation, tt is, in fact, a stationary solution of (5) i.e. can be obtained by solving 
the equation 

(8) 7zQ = 0. 

This homogeneous linear equation has non-trivial solutions because the matrix Q is singular 
by (6). An argument based on Perron-Frobenius theory for non-negative matrices ensures 
that TT is a unique strictly positive stochastic vector. Moreover, tt is stable so that for f —> oo 
indeed p{t) approaches tt i.e. we have limf^oo p{t) = n [68]. 

3.5. Estimation of Markov models from experimental data. Whereas the mathematical frame¬ 
work of aggregated Markov models was developed a short time after single channel data be¬ 
came available, the statistical estimation of these models is a topic of current research. Most 
commonly used are approaches based on Bayesian statistics. For a given time series Y of open 
and closed events recorded from an ion channel the conditional probability density f{Q\Y), 
known as the posterior density in the Bayesian framework, is used for determining a suitable 
Markov model with infinitesimal generator Q. Note that both Y and Q are considered as ran¬ 
dom variables, thus the posterior distribution quantifies how likely a model Q is under the 
condition that data Y have been observed. Direct calculation of the posterior /(Q|Y) is an¬ 
alytically intractable and computationally prohibitive but efficient approaches for maximum 
likelihood estimation (MLE) i.e. estimating 

(9) Q = argmaxQ/(Q|Y) 

were published in the 1990s [60, 61, 19]. Software implementations of these methods have 
been made available freely for academic use. Currently, the methods by Qin et al. [60, 61] can 
be obtained under the name QUB as standalone GUI applications at http: / / www. qub. buffalo. 
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Figure 3. Two examples for marginal distributions of rate constants, (a) shows 
a histogram with a well-defined mean ^ and a low standard deviation a which 
indicates a low level of parameter uncertainty whereas the histogram in (b) 
shows a complex multi-modal distribution which shows that multiple values 
of the rate constants are capable of representing the data. 


edu/. DCPROGS based on Colquhoun et al. [19] is still under active development and the source 
code of the most recent version has been published on github: https : //github. com/DCPROGS. 

An alternative approach to maximum likelihood estimation has been pursued since the late 
1990s. The aim of Markov chain Monte Carlo (MCMC) is to approximate the posterior den¬ 
sity f{Q\y) by sampling. MCMC enables us to randomly generate a sequence of 

models such that the expected frequency of a model within this sequence is as large as the 
density /(Q*^| Y). Thus, by generating a sufficient number of samples, the posterior /(Q|Y) is 
approximated. 

The early method by Ball et al. [4] for estimation of a Markov model Q depends on a suit¬ 
able idealisation of discretely sampled measurements to continuous open and closed times. 
This leads to a difficult statistical problem that has been discussed widely in the ion channel 
literature as the "missed events" problem. Rosales and colleagues were the first to propose 
a method that directly uses the discrete measurements and thus does not require further ide¬ 
alisation of the data [65, 64]. Their algorithm estimates a discrete-time Markov model which 
describes the transition probabilities between states during a sampling interval rather than 
the so-called infinitesimal generator Q. Gin et al. [36] were the first to propose a method for 
estimating Q from discretely-sampled data, their method was extended to models with ar¬ 
bitrary numbers of open and closed states by Siekmann et al. [74] and Siekmann et al. [70]. 
The current version of the software implementation of this method is available on github: 
https: //github. com/merlinthemagician/ahmm. git. For an overview of various approaches 
to statistical modelling based on single-channel data, see Gin et al. [38]. 

The crucial advantage of MGMG methods over MLE approaches is that uncertainties can be 
comprehensively understood by analysing the posterior /(Q|Y). Already marginal distribu¬ 
tions for individual rate constants (Figure 3) are helpful for localising and quantifying uncer¬ 
tainties within a model Q. But even more can be gained by analysing statistical relationships 
between combinations of model parameters as, for example, demonstrated by Siekmann et al. 
[70]. An important drawback of aggregated Markov models is non-identifiability i.e. model 
structures whose parameters cannot be inferred unambiguously from experimental data. Un¬ 
fortunately, non-identifiable aggregated Markov models have not been completely classified 
[32, 33, 11]. But non-identifiability can at least be detected by analysing the posterior distribu¬ 
tion /(Q|Y) [70]. Thus, MGMG allows us to disentangle different causes of model uncertainty 
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because it enables us to distinguish between parameter uncertainties due to insufficient or 
noisy data from pathologies in the structure of the model itself. 

3.6. The Ullah et al. model. A common approach for selecting a model structure for an ion 
channel model (which goes back at least to the classical model by De Young and Keizer [25]) 
is to identify the states of the Markov model with different chemical states of the channel 
protein. As explained in Section 3.1, the IP 3 R has various binding sites that allow specific 
ligands such as Ca^+ and IP 3 to bind to the channel protein and induce conformational changes 
of its three-dimensional structure. To account for this, model states are distinguished by how 
many particles of each ligand are bound to the channel. This assumption not only determines 
the state set of the model but also the possible transitions between states—in each state we can 
either bind a ligand to a free binding site or remove a ligand from an occupied binding site. 
The dynamics of binding and unbinding of ligands is modelled by the law of mass action so 
that, in principle, the model is completely specified by the number of binding sites for each 
ligand. However, in practice, such a model would be heavily overparameterised when fitted 
to experimental data, so it is necessary to simplify the model. 

To illustrate this with an example, consider the recent model by Ullah et al. [84] which is 
representative for this approach. The model states in Figure 2a are arranged in a grid whose 
horizontal axis shows how many Ca^^ molecules are bound to the channel and whose vertical 
axis indicates how many IP 3 binding sites are occupied. Thus, the position within the grid re¬ 
flects for a specific model state how many Ca^+ ions and how many IP 3 molecules, respectively, 
are bound to the channel. For example, neither Ca^+ nor IP 3 are bound to the state Cqq in the 
lower left corner whereas two Ca^+ and four IP 3 binding sites are occupied for the states C\^, 
O 24 , and O^. This is also indicated by the subscript indices—the first digit stands for the 
number of Ca^^ ions whereas the second digit accounts for the number of IP 3 molecules bound 
to the channel. Figure 2a shows that of the 20 possible combinations of occupying Ca^+, ATF 
and IF 3 binding sites only a subset of eight appears in the model. This considerable reduction 
is due to the removal of "low occupancy states"—Ullah et al. [83] developed a perturbation 
theory approach that allows them to omit states with low stationary probabilities while at the 
same time accounting for the delays caused by passing through these states. 

The model is constructed in an iterative four step process integrating several sources of data. 
In a first step, Ullah et al. [84] use Ca^+ and IF 3 dependency of the average open probability po 
in order to determine a minimal set of model states. By optimising an Akaike information 
criterion (AIC) score function, a model with five closed, Cqq, Cq 4 , C 24 , C 32 and C 34 , and one 
open state, O 24 , was selected as the best fit for the po data. 

In a second step, the ligand-dependent average probabilities tt^, and tz^ of being in 
modes characterised by three different levels of activity as well as the open probabilities in 
each mode (Pq, Pq and Pq) are used for assigning each of the six model states with a mode. At 
this step, some additional states are added because, for example, the state Cq 4 must exist both 
in the low (CQ 4 ) as well as the intermediate mode (CQ 4 ) in order to get a good fit to the data. To 
account appropriately for the Ca^+ dependency of Pq, the open probability in the intermediate 
mode, an additional state O 44 had to be introduced. 

In the first two steps, Ullah et al. [84] use stationary probabilities in order to determine which 
states should appear in the model without considering transitions between states. In step 3 the 
authors infer the transitions that are needed to account for the average sojourn times t^, 
and in the three modes whereas in step 4, data on the IF 3 R response to rapid changes 
in Ca^^ and IF 3 (latencies) is used for determining the remaining transitions. Two additional 
states, C 20 and C 3 Q are introduced in order to account for the latency data. 

Until this point, data is only used for determining the model structure but not for parameter 
estimation. The model is finally parameterised using the latency data from Mak et al. [54] or 
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a combination of these data and single-channel time series obtained at three different constant 
Ca^+ concentrations. 

3.7. Siekmann et al. "Park-Drive" model. The main aims of the modelling study by Siek- 
mann et al. [73] were first to account for switching between an inactive "park" and an active 
"drive" mode observed in the data set by Wagner and Yule [ 86 ]. As mentioned by Mak et al. 
[54] and Foskett and Mak [30], Mak and Foskett [53], the importance of modal gating is well- 
recognised and the implications for not appropriately capturing the timescale separation of 
fast opening and closing and slower switching between different activity levels is obviously 
unsatisfactory from a modelling point of view. 

Second, these data provided the possibility to build a model of two different mammalian 
isoforms of the IP 3 R, type I and type IIIP 3 R. In addition to a comparative study of type I and 
type IIIP 3 R, these data also include ligand-dependency of ATP in addition to IP 3 and Ca^+. 

Third, Siekmann et al. [73] followed a primarily statistical approach to inference, rather than 
deriving the model from a binding scheme as the model by Ullah et al. [84] discussed above. 
Based on the experience of the earlier study by Gin et al. [37] where similar data could be fitted 
satisfactorily by a model with four states and only one ligand-dependent pair of rate constants, 
the number of parameters required to account for binding of IP 3 , Ca^+ and ATP were likely to 
lead to a highly overparameterised model. 

Due to these considerations, Siekmann et al. [73] made the inactive "park" and the active 
"drive" mode the construction principle of their model. In a first step, Markov models repre¬ 
senting the stochastic dynamics for these two modes were constructed based on representative 
segments of the time series data that were characteristic for one of the two modes. Models 
with different numbers of states and model structures were fitted to these segments using the 
method by Siekmann et al. [74,70]. It was observed that the best fits for either of the two modes 
across all combinations of ligands available in the large data set by Wagner and Yule [ 86 ] were 
quantitatively similar. This strongly suggested (consistent with lonescu et al. [44]) that the 
dynamics within park and drive modes are ligand-independent and that ligand-dependent 
regulation of IP 3 R activity is achieved by varying the prevalence of park or drive mode. 

In a second step after both park and drive mode had been modelled separately, a model 
of the ligand-dependent switching between the ligand-independent modes was constructed. 
The structure for the full Park-Drive model (Figure 2b) was found by connecting the Markov 
models of park and drive mode obtained previously with a pair of transition rates. Due to 
the infrequent switching between park and drive mode observed in the data it was decided 
that adding more than a single pair of transition rates was statistically unwarranted. The full 
Park-Drive model was then fitted to time series for all combinations of ligands of the study 
by Wagner and Yule [ 86 ] . The results of these fits established the ligand-dependency of modal 
gating by the IP 3 -, Ca^^- and ATP-dependent variation of the two transition rates. 

Probably the most important result of this study is that only models that take into account 
modal gating are able to accurately capture IP 3 R kinetics. A channel whose kinetics is restricted 
to an active and an inactive mode can produce intermediate activity only by switching between 
both modes. Thus, a model that is not capable of switching between different levels of activity 
is misleading because it produces a constant open probability instead of alternating between 
highly different open probabilities. However, Cao et al. [13] showed that accounting for modal 
gating alone was insufficient for modelling stochastic Ca^^ release events (puffs) that arise 
from the interactions of a few IP 3 R channels. This study showed that the Park-Drive model has 
to be augmented by latency data [54] in order to account for the delayed response of individual 
channels to changes in ligand concentrations. 

Constructing the Park-Drive model based on the two modes proved very useful in the study 
by Cao et al. [14]. The authors iteratively reduced the Park-Drive model to a two-state model 
that only approximates the dynamics of opening and closing within the modes and focuses on 
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the level of activity determined by the relative prevalence of the modes. This further empha¬ 
sises that switching between park and drive mode rather than stochastic dynamics within the 
modes is the most important mechanism of IP 3 R regulation. 

3.8. Comparison of type I and type IIIP 3 R. The experimental study by Wagner and Yule [ 86 ] 
not only investigated the IP 3 R under a wide range of ligand conditions but also contrasted the 
behaviour of type I and type IIIP 3 R. In the models for type I and type IIIP 3 R constructed by 
Siekmann et al. [73] at a first glance the similarities between both subtypes are probably more 
obvious than the differences. First of all, it is striking that both IP 3 R subtypes can not only be 
represented in the same model structure but that active and inactive modes in both channels 
are nearly identical. This indicates that both subtypes have the same modes and that their 
differences are entirely due to differences in modal gating. 

One difference is that type IIIP 3 R responds more sensitively to IP 3 , in contrast to type IIP 3 R. 
The most important differences between both subtypes was found to be ATP regulation, see 
Wagner and Yule [ 86 ], Siekmann et al. [73] for details. 

4. Using data-driven IP 3 R models in calcium dynamics 

So far we have focused on the dynamics of individual IP 3 RS. In order to investigate the role 
of IP 3 RS in calcium dynamics we will now consider the interaction of IP 3 RS within a cluster. 

4.1. Modeling calcium puffs using the Park-Drive IP 3 R model. There is a large literature on 
stochastic models of calcium puffs for which we refer to the recent review by Rudiger [ 66 ]. 
Here we present a simple model based on the Park-Drive model [73] which is based on the 
following assumptions: 

• The ER contains sufficiently high [Ca^+]i to keep a nearly constant Ca^+ release rate 
through a cluster of IP 3 R [85]. Thus, ER [Ca^+]i dynamics is not explicitly modeled. 

• Ca^+ fluxes through the cell membrane have little effect on the very localised Ca^+ 
puffs far from cell membrane. 

• We compartmentalise our model to capture heterogeneity within a cluster of IP 3 RS. We 
assume that sufficiently far away from individual channels we have a homogeneous 
basal Ca^+ concentration c =[Ca^+]i that slowly responds to the total Ca^+ flux /ipjR 
through all IP 3 R channels. In the vicinity of an open IP 3 R channel this basal concentra¬ 
tion c is elevated by a constant C;,; once the channel closes it instantaneously equilibrates 
to the basal concentration c. 

Eurthermore, Ca^+ buffers are not considered except a Ca^^ fluorescence dye. With these 
assumptions, the model is given as follows, 

(10) ^ = /1P3R + /leak — ^ — kon{B — b)c -|- k^fib 

db 

(11) — = kon{B - b)c - koifb 

where V^c /(c -|- fC^) models the flux (mainly via diffusion and SERCA) removing Ca^+ from 
the puff site, /leak represents Ca^+ leak current from the ER for stabilising the resting [Ca^+ji 
of O.lfiM (a typical value). B and b represent the total dye buffer concentration and Ca^^- 
bound dye buffer concentration respectively, and the buffering process follows the mass action 
kinetics. /iPgR is the Ca^+ flux through open IP 3 R, which is modeled by the production of a 
constant release flux rate (Jr) and number of open IP 3 R channels (No), i.e. /iPgR = JrNg. Each 
open IP 3 R will equally contribute to the elevation of cluster [Ca^+ji, c. Note that the actual 
[Ca^+]i modulating each IP 3 R is either c (when it is in closed states) or c + Ch (when it is in open 
states). Parameters values are Jr = 200 fiM/s, V/ = 4000 jiM/s, = 12 }iM, /leak = 33 }iM/s, 
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B = 20 }iM, kon = 150 /Coff = 300 and = 120 jiM [13]. The cluster is assumed 

to contain 10 IP 3 R channels. 

The Park-Drive IP 3 R model is used to simulate IP 3 R state and coupled to the deterministic 
equations via a hybrid-Gillespie method [67]. However, the puff model based on the Park- 
Drive model fails to reproduce nonexponential interpuff interval (IPI) distribution due to the 
sole use of stationary single channel data (i.e. Ca^+ is fixed during measurement) in IP 3 R 
model construction. This does not allow the model to capture the transient single charmel 
behaviour when Ca^+ experiences a rapid change [54, 13]. Thus, the Park-Drive model is 
modified by incorporating time-dependent inter-mode transitions so that the transient single 
channel behaviour is captured. In detail, the transition rates (^24 arid q ^2 are changed from 
constants to functions of four newly-introduced gating variables, 

(12) q24: = fl24 + 1 ^ 4(1 — WJ24^24) 

(13) q42 = <^42 + V42WJ42^42 
where m 24 , h 24 , ^42 and /z 42 are gating variables obeying 

dG 

(14) — = Ag(Goo — G), G = WJ24,1224/42/^42• 

Goo is the steady state which is a function of channel-sensed Ca^^ and IP 3 concentrations and 
is determined by stationary single channel data (i.e. the Park-Drive model). Ac is the rate 
at which the steady state is approached. This is based on the fact that a IP 3 R channel cannot 
immediately reach its steady state upon a transient change in Ca^+ concentration [54]. The val¬ 
ues of Ag for 24224/1224 and 22242 are found to be large so that the three gating variables could be 
approximated by their steady states i.e. G = Goo/ a method called quasi-steady-state approx¬ 
imation. However, A;,^^ at low [Ca^+ji should be very small, reflecting a very slow recovery 
of IP 3 R from high Ca^+ inhibition [54]. Note that when A/j^^ is sufficiently large, quasi-steady- 
state approximation applies and the modified IP 3 R model reduces to the original Park-Drive 
model. Details about the functions and parameters can be seen in [13]. 



Figure 4. A simulation result of calcium puffs. F/Fq represents the ratio of b 
to its resting value. IP 3 concentration is O.ljiM. Adopted from [13]. 

An example of simulation results using the modified Park-Drive model is give in Fig. 4. 
The waiting time between two successive puffs (or interpuff interval, IPI) is a key statistics 
to quantify the underlying process governing the emergence of puffs. Fig. 5 shows that, as 
A;j ^2 low [Ca^+]i increases, the IPI distribution changes from nonexponential to exponential, 
demonstrating that the missing slow time scale in the original Park-Drive model is very crucial 
to explain the inhomogeneous Poisson process governing puff emergence found by (Thurley 

14 










et al. [81]). The IPI distributions were generated by fitting the probability density function pro¬ 
posed by Thurley et al. [81] to the simulated IPI histograms [13]. The proposed IPI distribution 
is 

(15) P = A(1 - 

where t represent IPI. A is the puff rate, a measure of the typical IPI (similar to average puff 
frequency), and ^ is the recovery rate. 



Figure 5. Dependence of IPI distribution on A;^^^ (indicated in the legend) at 
low [Ca^+]i. IP 3 concentration is O.lpM. Adopted from [13]. 

Hence, this example shows the particular importance of considering both stationary and 
nonstationary data when constructing an IP 3 R model. However, even if a model is constructed 
based on both data sets, it could also fail to reproduce Ca^+ puffs. One example is the Ullah 
model [84] as introduced in Section 3.6. A model simulation using the same puff model (10), 
(11) with the Ullah model is given in Fig. 6 where the Ca^+ signal behaves very irregularly and 
no puffs are clearly detected. 

4.2. The role of modal gating of IP 3 R in modulating calcium signals. The Park-Drive model 
(and its modified version) has the feature that IP 3 R exist in two different modes, each of which 
contains multiple states, some open, some closed. Intermode transitions are important for 
modulating Ca^+ signals because of their ligand- and time-dependent property. However, 
structure within each mode may also have substantial contribution to the formation of dif¬ 
ferent Ca^+ signals. Here, we examine the relative importance of intermode and intramode 
transitions using model reduction methods. By reducing the 6 -state IP 3 R model to a 2-state 
open/closed model, we will remove the intramodal structure, and a direct comparison be¬ 
tween the statistics generated by the two IP 3 R models will show the importance of intramodal 
structure. 

The model reduction takes the following steps: 

• The low probabilities of Ci , C 3 and O 5 (sum of which is less than 0.03 for any [Ca^+ji) 
means that the IP 3 R either rarely visit those states or have very short dwell time in 
those states. This allows to completely remove the three states from the 6 -state model. 

• Transitions q 26 and q (,2 are far larger (about two orders of magnitude) than and (^ 42 . 
By taking a quasi-steady state approximation to the transition between C 2 and Oe, we 
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Figure 6. A simulation result of calcium puffs using the Ullah IP 3 R model. IP 3 
concentration is O.ljiM. y-axis values indicate the number of IP 3 R channels in 
corresponding states. Parameter values for the puff model remain the same. 


have 0(, = € 2 ^ 26 /Combining C 2 and Oe to be a new state D, i.e. D = C 2 + Oe, the 
6-state model becomes a 2-state model, where D represents a partially open state with 
Ca^+ flux through the channel decreased by a factor of q 2 e/ {qei + ^26)- Moreover, ^^24 
needs to be rescaled by qsi/ [qei + due to the quasi-steady state approximation so 

that the effective closing rate is q2A%2/{qei + qie)- 
Fig. 7 shows the distributions of interpuff interval, puff duration and amplitude generated 
by using the 6 -state IP 3 R model (the Park-Drive model) and the reduced 2-state model. Re¬ 
ducing the intramodal structure does not qualitatively change the distributions but may lead 
to quantitative difference, which could be caused by missing open state O 5 that significantly 
contributes to the fluctuations of basal level of [Ca^+]i. However, if the IP 3 R channel is not 
very sensitive to small fluctuations of basal [Ca^+]i, the quantitative difference is significantly 
reduced [14]. Thus, the fundamental process governing the generation of Ca^+ puffs and oscil¬ 
lations is primarily controlled by the modal structure but not the intramodal structure which 
improves the model fitting to the single-channel data. 

5. Conclusions 

The IP 3 R plays a major role in CICR. For this reason, more and more aspects of its behaviour 
have been investigated by experiments. It usually turned out that new types of data had to be 
explicitly included in a model to account for them. For example, in early models such as the 
de Young-Keizer model [25], the rate constants were determined by fitting to the po observed 
at different calcium concentrations. But it soon became obvious that models parameterised 
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Figure 7. Comparison of interpuff interval, puff duration and amplitude be¬ 
tween the 6 -state IP 3 R model (the Park-Drive model) and the reduced 2-state 
model. 199 samples for each model were used to generate A and 200 samples 
for B and C. Interpuff interval distributions were fit by using Eq.15 proposed by 
Thurley et al. [81]. Puff amplitude distributions were fit by normal distribution. 


with po data could not be used for extrapolating the channel kinetics, i.e. the stochastic open¬ 
ing and closing. See Sneyd and Falcke [77] or Ullah et al. [84] for a more detailed explanation 
why it is impossible to infer kinetics from the ligand dependency of the open probability po- 

Just as kinetics cannot be inferred from po it turned out that the response of the IP 3 R to 
varying ligand concentrations cannot be predicted from data collected at constant ligand con¬ 
centrations. This was demonstrated by the next generation of models that were directly fitted 
to single-channel data, taking into account the stochastic process of opening and closing. The 
simplest assumption for integrating models for different ligand concentration is that the IP 3 R 
adjusts instantaneously. If this was true we could represent the channel kinetics appropriately 
by (to give a concrete example) simply replacing the model for the kinetics at 0.05 pM with 
the model for the kinetics at 0.2 pM calcium as we increase the calcium concentration. But 
Cao et al. [13] showed that only after taking into account rapid-perfusion data generated by 
Mak et al. [54] was the model of Siekmann et al. [73] capable of generating the correct puff 
distribution. 

It is important to note that taking into account more data does not necessarily have to lead 
to more complicated models. Instead, after taking into account that the simpler kinetics of 
modal gating should capture the part of the channel dynamics that is most important for the 
functional role of the IP 3 R in CICR, Cao et al. [14] were able to reduce the six-state model by 
Siekmann et al. [73] to a two-state model. Thus, after interpreting experimental data in the 
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right way, we are able to build models for the functional role of IP 3 R that are nearly as simple 
as the early phenomenological models. 

6. Future work 

After reviewing the current state of data-driven approaches to investigating the IP 3 R we 
would like to take a look at promising future directions. In order to address the particular 
importance of modal gating, Siekmann et al. [71] develop a novel hierarchical model structure 
that enables us to combine Markov models that represent the stochastic switching between 
modes with models that account for the characteristic opening and closing within different 
modes. Thus, models for both processes can be fitted separately (e.g. using the method by 
Siekmann et al. [74,70]) after analysing the data with statistical method presented by Siekmann 
et al. [72]. This allows us to build models for modal gating following a completely data-driven 
approach. 

More generally, we have compared two current models as representative examples for dif¬ 
ferent modelling approaches, the Ullah et al. [84] and the Park-Drive model [73, 13, 14]. Al¬ 
though both approaches ultimately meet in the middle, their different construction principles 
impose different requirements for future progress. From a statistical point of view, representa¬ 
tion of ligand interactions with a channel by mass action kinetics as in Ullah et al. [84] defines 
a sufficiently large search space of models. It is crucial to select from this search space an 
appropriately simplified model that is obtained by removing states of the full model in a con¬ 
sistent way. A method for model reduction is provided by Ullah et al. [83] and Ullah et al. [84] 
demonstrate how data can be used to statistically select from all possible simplified models. A 
central principle of the biophysical approach is to design models in a way that closely follows 
physical principles. In this context, the bond-graph approach to modelling ion channels by 
Gawthrop and Crampin [35], Gawthrop et al. [34] is highly relevant because it ensures that 
physical principles are enforced when choosing a model structure. 

For models that primarily focus on a statistically satisfying representation in a first instance, 
the model selection problem arises again but in the other direction. Rather than starting from 
a model structure determined by an underlying mass action model. Gin et al. [37] and Siek¬ 
mann et al. [73] iteratively increased the number of states in their model structure until fur¬ 
ther increasing the number of parameters appears statistically unwarranted. This process is 
time-consuming and may be computationally prohibitive if models exceed a certain number 
of states. Developing a method that is able to automatically compare models with an increas¬ 
ing number of states has proven to be difficult indicated by the few number of studies that 
have appeared on this subject after an early article on comparison of a finite number of models 
[41]. A promising new direction is the non-parametric Bayesian method developed by Hines 
et al. [39] which allows the authors to estimate the number of states within an ion-charmel data 
set. Determining the required number of open and closed states in a first step may increase 
efficiency because it restricts the class of models which have to be compared in a second step. 
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